Analysis and Making plot for paper

光合成実験

• 小型回流水槽に暗幕をし、明条件1h, 暗条件2hとした。
• 密度条件ごとに個体を変えた（流速条件では同じ）
• 実験開始前に、海水のDOを窒素ガスでばっ気することで、90%近くまで下げた。
• 光合成実験の際には、海水の空気中の酸素ガスの交換が発生する。そのため、ガス交換速度を流速条件ごとに算出し、補正した。

流速測定実験

• 流速はドップラー流速計（サンプリング間隔：200Hz、測定時間：3分）を用いて測定した。
• 測定位置は、藻場の直前、leading edgeから15 cm, 45 cm, 藻場直後とし、測定深度は、底面から6 cm, 藻場上部の境界（流速、密度によって異なる）、16 cmとした。

キャノピーの形状変化の測定

• キャノピーの流速による形状変化は、キャノピーの横方向からの投影面積で表現した。
• カメラで流速ごとに撮影し、imagejを用いて面積を算出した。

光量の測定

• 光量は約800\(\mu mol \: photons \: m^{-2}\:s^{-1}\)の一定。
• 藻場底面の流速、密度ごとの光量は、ペンダントロガーにて測定（lux）

溶存酸素濃度の生データの確認

それぞれの条件ごとに、測定初期の乱れを除去するため900秒を捨てる。変動速度測定に使用するのは、開始位置から40分間とする。 TRUEが明条件、FALSEが暗条件を示す。

## # A tibble: 1,968 x 11
## # Groups:   Hz, Light, Density [48]
##       ID    Hz Light Elapsed_time Density    DO  Temp Barometer   Sal
##    <dbl> <dbl> <lgl>        <dbl> <fct>   <dbl> <dbl>     <dbl> <dbl>
##  1     3   0.5 FALSE            0 High     6.82  21.3      101.  34.6
##  2     3   0.5 FALSE            1 High     6.81  21.3      101.  34.6
##  3     3   0.5 FALSE            2 High     6.82  21.3      101.  34.6
##  4     3   0.5 FALSE            3 High     6.81  21.3      101.  34.6
##  5     3   0.5 FALSE            4 High     6.82  21.3      101.  34.6
##  6     3   0.5 FALSE            5 High     6.8   21.4      101.  34.6
##  7     3   0.5 FALSE            6 High     6.81  21.3      101.  34.7
##  8     3   0.5 FALSE            7 High     6.81  21.3      101.  34.6
##  9     3   0.5 FALSE            8 High     6.81  21.4      101.  34.6
## 10     3   0.5 FALSE            9 High     6.81  21.4      101.  34.6
## # ... with 1,958 more rows, and 2 more variables: WW_sum <dbl>,
## #   Seaweed_volume_m3 <dbl>

#溶存酸素濃度の変化速度を算出するための関数作成
Glmmodel = function(y, x){
  summary(glm(y ~ x, family = Gamma(link = "log")))
}　　　　　　　　　　　　　　　　　　　　　　　　

Slope = function(x){　　　　　　　　　　　　　　　
  coef(x)[2]
} 

Slope_error = function(x){　　　　　　　　　　　　　　　
  coef(x)[4]
} 

#溶存酸素変化速度を算出算出
DO_analysis = Trident %>%
  group_by(Hz, Light, WW_sum, Seaweed_volume_m3, Density) %>%
  nest() %>% 
  mutate(Model = map(data, ~ glm(DO ~ Elapsed_time, data = ., family = Gamma(link = "log")))) %>% 
  mutate(Summary = map(Model, ~summary(.))) %>% 
  mutate(Slope = map_dbl(Model, ~Slope(.)),
         Stderr = map_dbl(Model, ~Slope_error(summary(.))))

#湿重量あたりの溶存酸素変化速度を算出
DO_analysis_WW = DO_analysis %>% 
　mutate(Slope = Slope / WW_sum, Stderr = Stderr / WW_sum) 

## # A tibble: 48 x 10
##       Hz Light WW_sum Seaweed_volume_… Density data  Model Summary    Slope
##    <dbl> <lgl>  <dbl>            <dbl> <fct>   <lis> <lis> <list>     <dbl>
##  1   0.5 FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -4.15e-7
##  2   0.5 TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  2.32e-6
##  3   1   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -5.82e-7
##  4   1   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  2.55e-6
##  5   2   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -4.68e-7
##  6   2   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  2.96e-6
##  7   4   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -5.04e-7
##  8   4   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  3.26e-6
##  9   6   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -2.50e-7
## 10   6   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  3.12e-6
## # ... with 38 more rows, and 1 more variable: Stderr <dbl>

溶存酸素変化速度の補正：海水と空気中のガス交換速度にて

　小型回流水槽で実験する際、測器の取り付けや実験の効率を考えて、観測部分水路の上部は開放した状態で実験を行った。そのため、流速が速くなるにつれて、海水-空気間の酸素ガスの交換が活発となり、溶存酸素濃度の変動に影響を与えると考えられる。そこで、海藻を入れていない状態で、窒素ばっ気により、溶存酸素濃度を低下させた海水の溶存酸素濃度の経時変化を測定し、流速条件ごとのガス交換速度を算出した。
　実験を行った流速条件は6, 0.5, 4, 8, 20 Hzとした。

補正に使用するデータの確認。

#大気-海水間の酸素フラックスの算出-------------------------------------------------------------
#光合成実験3で測定したデータを使用する。回流水槽内に海藻は入っていない。
Trident_flux = read_csv("../回流水槽光合成実験3/Modified_data/DO_air_water_Flux_data.csv") 
print(Trident_flux, n = 10)
ggplot(Trident_flux)+
  geom_point(aes(x = Elapsed_time, y = DO, color = "DO"))+
  geom_point(aes(x = Elapsed_time, y = DO_saturation, color = "DO_saturation"))+
  facet_wrap("Hz")

## # A tibble: 3,177 x 11
##    Time                ID    Barometer    DO  Temp   Sal
##    <dttm>              <chr>     <dbl> <dbl> <dbl> <dbl>
##  1 2018-05-06 06:15:27 calb…      101.  6.88  18.9  34.7
##  2 2018-05-06 06:16:27 calb…      101.  6.88  19    34.7
##  3 2018-05-06 06:17:27 calb…      101.  6.88  19    34.7
##  4 2018-05-06 06:18:27 calb…      101.  6.89  19    34.7
##  5 2018-05-06 06:19:27 calb…      101.  6.89  19    34.7
##  6 2018-05-06 06:20:27 calb…      101.  6.89  19    34.7
##  7 2018-05-06 06:21:27 calb…      101.  6.9   19    34.7
##  8 2018-05-06 06:22:27 calb…      101.  6.9   19    34.7
##  9 2018-05-06 06:23:27 calb…      101.  6.9   19    34.7
## 10 2018-05-06 06:24:27 calb…      101.  6.91  19    34.7
## # ... with 3,167 more rows, and 5 more variables: Date <dttm>, Hz <dbl>,
## #   DO_saturation <dbl>, DO_diff <dbl>, Elapsed_time <dbl>

　以上のデータを以下の式に当てはめる。 　

　\(DO_{t}\:=\:C\:*\:exp(-k\:*\:t)\:+\:DO_{sat}\)
　\(DO_{t}\):経過時間tにおける溶存酸素濃度
　\(C\):t = 0における溶存酸素濃度と飽和溶存酸素濃度の差
　\(k\):溶存酸素濃度のガス交換による変化速度
　\(DO_{sat}\):経過時間tにおける飽和溶存酸素濃度
　

#データを用いてKの算出（ｋが負になると酸素の放出、正は吸収）-----------------------------
library(nlstools)
model = function (k , C, t, DO_saturation){
  C * exp(-k * t) + DO_saturation
}

k = Trident_flux %>%
  filter(!ID == "ECSER09") %>% 
  group_by(Hz) %>% 
  nest() %>% 
  mutate(nlsout = map(data, function(x){nls(DO~model(k0,C0,Elapsed_time, DO_saturation), data = x, start = list(k0=0.11, C0=-4))}),
         k = map_dbl(nlsout, function(x){coef(x)[[1]]}))

#当てはまり確認
k
k$nlsout

## # A tibble: 5 x 4
##      Hz data                nlsout           k
##   <dbl> <list>              <list>       <dbl>
## 1  20   <tibble [241 × 10]> <S3: nls> 0.00554 
## 2   8   <tibble [395 × 10]> <S3: nls> 0.00213 
## 3   4   <tibble [665 × 10]> <S3: nls> 0.000864
## 4   0.5 <tibble [318 × 10]> <S3: nls> 0.000138
## 5   6   <tibble [620 × 10]> <S3: nls> 0.00244 
## [[1]]
## Nonlinear regression model
##   model: DO ~ model(k0, C0, Elapsed_time, DO_saturation)
##    data: x
##        k0        C0 
##  0.005539 -0.655192 
##  residual sum-of-squares: 0.006818
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 0.00000138
## 
## [[2]]
## Nonlinear regression model
##   model: DO ~ model(k0, C0, Elapsed_time, DO_saturation)
##    data: x
##        k0        C0 
##  0.002126 -0.732538 
##  residual sum-of-squares: 0.01068
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 0.000008053
## 
## [[3]]
## Nonlinear regression model
##   model: DO ~ model(k0, C0, Elapsed_time, DO_saturation)
##    data: x
##         k0         C0 
##  0.0008645 -0.7389470 
##  residual sum-of-squares: 0.01915
## 
## Number of iterations to convergence: 8 
## Achieved convergence tolerance: 0.00000008483
## 
## [[4]]
## Nonlinear regression model
##   model: DO ~ model(k0, C0, Elapsed_time, DO_saturation)
##    data: x
##         k0         C0 
##  0.0001376 -2.4797074 
##  residual sum-of-squares: 0.009746
## 
## Number of iterations to convergence: 7 
## Achieved convergence tolerance: 0.000009302
## 
## [[5]]
## Nonlinear regression model
##   model: DO ~ model(k0, C0, Elapsed_time, DO_saturation)
##    data: x
##        k0        C0 
##  0.002439 -4.372584 
##  residual sum-of-squares: 1.053
## 
## Number of iterations to convergence: 6 
## Achieved convergence tolerance: 0.00000157

すべてのHzのおいてｋを算出していないため、Hzとｋは直線回帰で近似できると仮定して、ｋを補足する

#算出した酸素フラックス速度定数kを持ちいて、実験中の酸素フラックスを見積もる:使用するモデルはパターン2のDO_saturation を変動モデル-----------------------------------------------------------------------------
#現在（2018-03-21現在）ではすべてのHzのおいてｋを算出していないため、Hzとｋは直線回帰で近似できると仮定して、ｋを補足する

Hz = data.frame(Hz = c(0.5, 1, 2, 4, 6, 8, 10, 15, 20))
k = k %>%
  full_join(Hz) %>% 
  arrange(Hz)

#Hz6の推定が悪いので、NAにする
k = k %>% mutate(k = ifelse(Hz == 6, NA, k))

#一次式で近似
k_pre = data_frame(Hz = k$Hz, k = predict(lm(k ~ Hz, k), data = data_frame(Hz = k$Hz),newdata=data.frame(Hz = k$Hz)))

ggplot()+
  geom_point(aes(x = Hz, y = k), data = k)+
  geom_line(aes(x = Hz, y = k), data = k_pre)

　kを用いて溶存酸素濃度を補正する

#kを用いてDO濃度を補正
DO_offset_Flux = Trident %>% 
  mutate(DO_saturation = gas_O2sat(S = Sal, t = Temp),
         DO_diff = DO_saturation - DO) %>% 
  left_join(k_pre, by = "Hz") %>%
  mutate(DO_Flux_min = DO_diff * k) %>% 
  group_by(ID, Hz, Light, Density) %>% 
  mutate(DO_Flux_min_cum = cumsum(DO_Flux_min),
         DO_offset_Flux = DO - DO_Flux_min_cum) %>%
  return()

print(DO_offset_Flux, n = 10)

ggplot()+
  # geom_point(aes(x = Elapsed_time, y = DO_Flux_min, color = "Flux"), data = DO_offset_Flux)+
  geom_point(aes(x = Elapsed_time, y = DO_offset_Flux, color = "補正後"), data = DO_offset_Flux)+
  # geom_smooth(aes(x = Elapsed_time, y = DO_offset_Flux, color = "offset_line"), method = "glm", data = DO_offset_Flux)+
  # geom_point(aes(x = Elapsed_time, y = DO_saturation, color = "Saturation"), data = DO_offset_Flux)+
  geom_point(aes(x = Elapsed_time, y = DO, color = "補正前"), data = DO_offset_Flux)+
  # geom_smooth(aes(x = Elapsed_time, y = DO, color = "raw_line"), method = "glm", data = DO_offset_Flux)+
  # geom_point(aes(x = Elapsed_time, y = DO_diff, color = "diff"), data = DO_offset_Flux)+
  # geom_point(aes(x = Elapsed_time, y = DO_Flux_min_cum, color = "cum"), data = DO_offset_Flux)+
  # geom_line(aes(x = Elapsed_time, y = k, color = "k"), data = DO_offset_Flux)+
  facet_grid(Density ~ Light + Hz)

## # A tibble: 1,968 x 17
## # Groups:   ID, Hz, Light, Density [48]
##       ID    Hz Light Elapsed_time Density    DO  Temp Barometer   Sal
##    <dbl> <dbl> <lgl>        <dbl> <fct>   <dbl> <dbl>     <dbl> <dbl>
##  1     3   0.5 FALSE            0 High     6.82  21.3      101.  34.6
##  2     3   0.5 FALSE            1 High     6.81  21.3      101.  34.6
##  3     3   0.5 FALSE            2 High     6.82  21.3      101.  34.6
##  4     3   0.5 FALSE            3 High     6.81  21.3      101.  34.6
##  5     3   0.5 FALSE            4 High     6.82  21.3      101.  34.6
##  6     3   0.5 FALSE            5 High     6.8   21.4      101.  34.6
##  7     3   0.5 FALSE            6 High     6.81  21.3      101.  34.7
##  8     3   0.5 FALSE            7 High     6.81  21.3      101.  34.6
##  9     3   0.5 FALSE            8 High     6.81  21.4      101.  34.6
## 10     3   0.5 FALSE            9 High     6.81  21.4      101.  34.6
## # ... with 1,958 more rows, and 8 more variables: WW_sum <dbl>,
## #   Seaweed_volume_m3 <dbl>, DO_saturation <dbl>, DO_diff <dbl>, k <dbl>,
## #   DO_Flux_min <dbl>, DO_Flux_min_cum <dbl>, DO_offset_Flux <dbl>

補正した溶存酸素濃度を用いて、変化速度を再計算

DO_analysis_k = DO_offset_Flux %>%
  group_by(Hz, Light, WW_sum, Seaweed_volume_m3, Density) %>%
  nest() %>% 
  mutate(Model = map(data, ~ glm(DO_offset_Flux~Elapsed_time, data = ., family = Gamma(link = "log")))) %>% 
  mutate(Summary = map(Model, ~summary(.))) %>% 
  mutate(Slope = map_dbl(Model, ~Slope(.)),
         Stderr = map_dbl(Model, ~Slope_error(summary(.))),
         Temp_mean = map_dbl(data, function(x){
           x %>%
             summarise(Temp = mean(Temp)) %>% 
             pull(Temp)}),
         Temp_sd = map_dbl(data, function(x){
           x %>%
             summarise(Temp = sd(Temp)) %>% 
             pull(Temp)}))

#湿重量あたりの溶存酸素変化速度を算出
DO_analysis_k_WW = DO_analysis_k %>% 
　mutate(Slope = Slope / WW_sum, Stderr = Stderr / WW_sum) 

補正後の溶存酸素変化量の確認。
呼吸速度は、0以下の値を取るように補正された。光合成速度は山形の変動を示した。

## # A tibble: 48 x 12
##       Hz Light WW_sum Seaweed_volume_… Density data  Model Summary    Slope
##    <dbl> <lgl>  <dbl>            <dbl> <fct>   <lis> <lis> <list>     <dbl>
##  1   0.5 FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -4.17e-7
##  2   0.5 TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  2.32e-6
##  3   1   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -6.09e-7
##  4   1   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  2.52e-6
##  5   2   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -5.47e-7
##  6   2   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  2.88e-6
##  7   4   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -6.37e-7
##  8   4   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  3.20e-6
##  9   6   FALSE   453.          0.00063 High    <tib… <S3:… <S3: s… -4.34e-7
## 10   6   TRUE    453.          0.00063 High    <tib… <S3:… <S3: s…  3.07e-6
## # ... with 38 more rows, and 3 more variables: Stderr <dbl>,
## #   Temp_mean <dbl>, Temp_sd <dbl>

溶存酸素変化速度の補正：水温による補正

補正は以下のアレニウスの式を用いて行う。Durocher and Allen (2012)を参照。

\(r^{'}\:=\:r\:*\:exp(-Ea\:/\:k\:(\:1/K\:-\:1/K0\:)\)

　\(r^{'}\):補正後（水温がk0のとき）の値
　\(r\):補正前（それぞれの値の水温）の値
　\(k\):ボルツマン定数
　\(Ea\):活性化エネルギー（海藻の状態や実験条件によって異なると考えられる）

海藻の生理的な応答がアレニウスの式の関係に従うと仮定し、K0のときの水温に合わせて補正する。 まず、補正前のデータもちいてEaを推定する。補正する水温は20 °C とした。

Eaの値が、以上に大きいため、この値で補正するには不適切であると考えられる。

bolzmann = 8.617e-5 #ボルツマン定数
K0 = 20 +  273.15 #補正水温K0は全ての条件で統一する。

#光条件ごとにEaを推定し、推定したEaでSlopeを補正
DO_analysis_k_WW_Ea = DO_analysis_k_WW %>% 
  group_by(Light) %>%
  mutate(K0 = K0,
         K = Temp_mean + 273.15,
         invK = (1 / K0 * bolzmann - 1 / K * bolzmann), 
         Slope = ifelse(Light == T, Slope, Slope*(-1)),
         Slslope = scale(log(Slope), scale = F)) %>% 
  nest() %>% 
  mutate(Model = map(data, function(x){lm(Slslope ~ invK, data = x) %>% summary()}),
         Ea = map_dbl(Model, function(x){coef(x)[2]}),
         data_temp = map2(data, Ea, function(x, y){
           x %>%
             mutate(Slope = Slope * exp(-y * invK))})) 

DO_analysis_k_WW_Ea
DO_analysis_k_WW_Ea %>% 
  select(Light, data_temp) %>% 
  unnest() %>% 
  ggplot()+
  geom_point(aes(x = Hz, y = Slope, shape = "補正前", color = Temp_mean),
             size = 3,
             data = DO_analysis_k_WW %>% 
               mutate(Slope = ifelse(Light == T, Slope, Slope*(-1))))+
  geom_point(aes(x = Hz, y = Slope, shape = "補正後", color = Temp_mean),
             size = 3)+
  facet_grid(Light ~ Density)

#Light,Density 条件ごとにEaを推定した場合。
DO_analysis_k_WW_Ea = DO_analysis_k_WW %>% 
  group_by(Light, Density) %>%
  mutate(K0 = K0,
         K = Temp_mean + 273.15,
         invK = (1 / K0 * bolzmann - 1 / K * bolzmann), 
         Slope = ifelse(Light == T, Slope, Slope*(-1)),
         Slslope = scale(log(Slope), scale = F)) %>% 
  nest() %>% 
  mutate(Model = map(data, function(x){lm(Slslope ~ invK, data = x) %>% summary()}),
         Ea = map_dbl(Model, function(x){coef(x)[2]}),
         data_temp = map2(data, Ea, function(x, y){
           x %>%
             mutate(Slope = Slope * exp(-y * invK))})) 

DO_analysis_k_WW_Ea
DO_analysis_k_WW_Ea %>% 
  select(Light, Density, data_temp) %>% 
  unnest() %>% 
  ggplot()+
  geom_point(aes(x = Hz, y = Slope, shape = "補正前", color = Temp_mean),
             size = 3,
             data = DO_analysis_k_WW %>% 
               mutate(Slope = ifelse(Light == T, Slope, Slope*(-1))))+
  geom_point(aes(x = Hz, y = Slope, shape = "補正後", color = Temp_mean),
             size = 3)+
  facet_grid(Light ~ Density)

## # A tibble: 2 x 5
##   Light data               Model                     Ea data_temp         
##   <lgl> <list>             <list>                 <dbl> <list>            
## 1 FALSE <tibble [24 × 15]> <S3: summary.lm> -286625485. <tibble [24 × 15]>
## 2 TRUE  <tibble [24 × 15]> <S3: summary.lm>  219348551. <tibble [24 × 15]>
## # A tibble: 6 x 6
##   Light Density data            Model                   Ea data_temp      
##   <lgl> <fct>   <list>          <list>               <dbl> <list>         
## 1 FALSE High    <tibble [8 × 1… <S3: summary.l…    -2.21e8 <tibble [8 × 1…
## 2 TRUE  High    <tibble [8 × 1… <S3: summary.l…     2.37e8 <tibble [8 × 1…
## 3 FALSE Low     <tibble [8 × 1… <S3: summary.l…    -1.85e7 <tibble [8 × 1…
## 4 TRUE  Low     <tibble [8 × 1… <S3: summary.l…    -3.46e8 <tibble [8 × 1…
## 5 FALSE Middle  <tibble [8 × 1… <S3: summary.l…    -1.74e8 <tibble [8 × 1…
## 6 TRUE  Middle  <tibble [8 × 1… <S3: summary.l…     2.23e8 <tibble [8 × 1…

次にEaを0.65で固定して補正する。
補正前後で大きな変化はない。水温の変動から考えても妥当か。

#全ての条件でEaを固定する。
FixedEa = 0.65 #1000にしてもほどんど変化なし。

DO_analysis_k_WW_FixedEa = DO_analysis_k_WW %>% 
  group_by(Light) %>%
  mutate(K0 = K0,
         K = Temp_mean + 273.15,
         invK = (1 / K0 * bolzmann - 1 / K * bolzmann), 
         Slope = ifelse(Light == T, Slope, Slope*(-1))) %>% 
  nest() %>% 
  mutate(Ea = FixedEa,
         data_temp = map2(data, Ea, function(x, y){
           x %>%
             mutate(Slope = Slope * exp(-y * invK))})) 

DO_analysis_k_WW_FixedEa
DO_analysis_k_WW_FixedEa %>% 
  select(Light, data_temp) %>% 
  unnest() %>% 
  ggplot()+
  geom_point(aes(x = Hz, y = Slope, shape = "補正前", color = Temp_mean),
             size = 3,
             data = DO_analysis_k_WW %>% 
               mutate(Slope = ifelse(Light == T, Slope, Slope*(-1))))+
  geom_point(aes(x = Hz, y = Slope, shape = "補正後", color = Temp_mean),
             size = 3)+
  facet_grid(Light ~ Density)

## # A tibble: 2 x 4
##   Light data                  Ea data_temp         
##   <lgl> <list>             <dbl> <list>            
## 1 FALSE <tibble [24 × 14]>  0.65 <tibble [24 × 14]>
## 2 TRUE  <tibble [24 × 14]>  0.65 <tibble [24 × 14]>

水温補正は、実験間の水温のばらつきが最大で約1°C 程であることから、それほど必要ないかもしれない。推定したEaの絶対値が大きくなったのは、Eaを算出する際に、密度条件や個体差によるSlopeの違いが、推定に反映されたからではないか。
よって水温補正は、固定したEaを用いた補正を使用する。

総光合成速度の算出

総光合成速度は、純光合成速度と総光合成速度の和で算出される。

NP = DO_analysis_k_WW_FixedEa %>% 
  select(Light, data_temp) %>% 
  unnest() %>% 
  filter(Light == TRUE) %>% 
  mutate(type = "NP") %>% 
  select(type, Hz, Density, WW_sum, Temp_mean, Seaweed_volume_m3, Slope)

RP = DO_analysis_k_WW_FixedEa %>% 
  select(Light, data_temp) %>% 
  unnest() %>% 
  filter(Light == FALSE) %>% 
  mutate(type = "RP") %>% 
  select(type, Hz, Density, WW_sum, Temp_mean, Seaweed_volume_m3, Slope)

GP = NP %>% 
  full_join(RP, by = c("Hz", "Density", "WW_sum", "Seaweed_volume_m3")) %>% 
  mutate(Slope = Slope.x + Slope.y, 
         type = "GP") %>% 
  select(type, Hz, Density, WW_sum, Seaweed_volume_m3, Slope)

#データの結合
V = 200 #L 回流水槽の容量
Slope_DO = bind_rows(NP, RP, GP) %>% 
  mutate(Slope = Slope*V, #単位をmg/l/min/gww -> mg/min/gww
         Slope = Slope*1000) #単位をmg/min/gww -> μg/min/gww

head(Slope_DO)

## # A tibble: 6 x 7
##   type     Hz Density WW_sum Temp_mean Seaweed_volume_m3 Slope
##   <chr> <dbl> <fct>    <dbl>     <dbl>             <dbl> <dbl>
## 1 NP      0.5 High      453.      21.7           0.00063 0.464
## 2 NP      1   High      453.      22.0           0.00063 0.505
## 3 NP      2   High      453.      21.6           0.00063 0.575
## 4 NP      4   High      453.      22.3           0.00063 0.640
## 5 NP      6   High      453.      21.9           0.00063 0.614
## 6 NP      8   High      453.      22.2           0.00063 0.588

流速データの読み込みとノイズの除去

vector_names = dir("../回流水槽光合成実験5/Data/Velocity/", pattern = "dat", full = T)

cnames_vectrino=c("time", "counter", "status",
                  "xvel", "yvel", "zvel", "z2vel",
                  "amp1", "amp2", "amp3", "amp4",
                  "snr1", "snr2", "snr3", "snr4",
                  "corr1", "corr2", "corr3", "corr4")

vector = lapply(vector_names, read_table , col_names=cnames_vectrino)

vector_test1 = vector("list", length(vector_names))
vector_test2 = vector("list", length(vector_names))
vector_test_result = vector("list", length(vector_names))

for(i in 1:length(vector_names)){
  # Try on a small subset of data
  # Trial run on the first 10000 data points
  vector_test1[[i]] = vector[[i]] %>% slice(1000:11000) #testのため一部抽出
  
  # Remove spikes
  # First clean up the data based on the SNR and CORR
  x = names(vector_test1[[i]] %>% select(dplyr::contains("vel")))
  vector_test1[[i]][,x]=snr_corr_filter(vector_test1[[i]])
  
  # Next run nk_despike() to clean it all up.
  vector_test2[[i]] = vector_test1[[i]] %>% mutate(tau = 1:length(xvel),
                                                   u = zoo::na.locf(xvel,na.rm=F), #NAを直近のNA以外の数値に置き換える
                                                   v = zoo::na.locf(yvel,na.rm=F),
                                                   w = zoo::na.locf(zvel,na.rm=F))
  
  vector_test_result[[i]] = vector_test2[[i]] %>% mutate_at(vars(u,v,w), nk_despike)
}

names(vector_test_result) = paste("", str_extract(vector_names, "[0-9]{1,2}_[0-9]{1,2}_[0-9]{1,3}_[HLM]"), sep = "")

vector_test_result2 = vector_test_result %>%
  bind_rows(.id = "location") %>%
  mutate(location2 = location) %>%
  tidyr::separate(location2, into=c("Hz", "Depth", "Distance", "Density"), sep = '_') %>%
  select(Hz, Depth,
         Distance, Density,
         time,
         xvel, yvel, zvel, z2vel,
         tau, u, v, w) %>%
  mutate(Depth = as.numeric(Depth),
         Hz = as.numeric(Hz)) %>% 
  mutate(Distance = as.numeric(Distance)) %>% 
  arrange(Hz, Distance, Depth, Density) %>%
  tbl_df()

vector_all = vector_test_result2 %>%
  mutate(Hz = ifelse(Hz == 5, 0.5, Hz)) %>% 
  mutate(Density = recode(Density, H = "High", M = "Middle", L = "Low"))

#保存
write.csv(vector_all, "../回流水槽光合成実験5/Modified_Data/vector_all.csv")

流動の統計量の算出

　算出流動に関する統計量は以下とする。

• 実験条件、測定値ごとの平均値: \(\bar{u}\)
  
• 残差: \(u^{\prime}=u-\bar{u}\)
  
• 残差の二乗平均: \(\bar{u}^{\prime}=\sum(u^{\prime})^{2}\:/\:n\)
  
• TKE(乱流エネルギー): \(TKE=0.5*(\bar{u}^{\prime2}+\bar{v}^{\prime2}+\bar{w}^{\prime2})\)
  
• Reynold stress: \(\bar{RS}= -\sum u^{\prime}*w^{\prime}\:/n\)
  　

vector_all$Density = factor(vector_all$Density, levels = c("High", "Middle", "Low"))　#密度の順位変更

vector_mean = vector_all %>%
  group_by(Hz, Depth, Distance, Density) %>%
  summarise_at(vars(xvel, yvel, zvel, u, v, w),funs(mean, sd)) 

vector_zansa = vector_all %>%
  left_join(vector_mean, by=c("Depth", "Distance", "Hz", "Density")) %>%
  mutate(zansaX = u - u_mean, zansaY = v - v_mean, zansaZ = w - w_mean)

vector_zansa_mean = vector_zansa %>%
  group_by(Depth, Distance, Hz, Density) %>%
  summarise(zansaX_mean = mean(zansaX^2),
            zansaY_mean = mean(zansaY^2),
            zansaZ_mean = mean(zansaZ^2),
            stress = mean(zansaX * zansaZ) * (-1),
            TKE = 0.5*(zansaX_mean + zansaY_mean + zansaZ_mean))

#藻場直前の測定値の鉛直平均を上流流速とする。
upstream_u_mean = vector_mean %>%
  filter(Distance == 26) %>%
  select(Hz, Depth, Density, upstream_u_mean = u_mean) %>% 
  group_by(Hz, Density) %>% 
  summarise(upstream_u_mean = mean(upstream_u_mean))

vector_mean = vector_mean %>% left_join(upstream_u_mean, by = c("Hz", "Density"))
vector_zansa_mean = vector_zansa_mean %>% left_join(upstream_u_mean, by = c("Hz", "Density"))

#測定深度を藻場の中、藻場の境界、藻場の上にグループ分けする。
vector_mean2 = vector_mean %>%
  mutate(Depth2 = ifelse(Depth >= 7 & Depth <= 13, "edge", as.character(Depth)),
         Depth2 = recode(Depth2, "16" = "above", "6" = "in")) 

vector_zansa_mean2 = vector_zansa_mean %>%
  mutate(Depth2 = ifelse(Depth >= 7 & Depth <= 13, "edge", as.character(Depth)),
         Depth2 = recode(Depth2, "16" = "above", "6" = "in")) 

光データの編集

光（lux）は、ペンダントロガーを２つを、 leading edgeから15 cm, 45 cmの底面に設置して測定した。インターバルは１秒。

#データの読み込み---------------------------------------------------------------------------------------------
fnames = dir("../回流水槽光合成実験5/Data/Light/", pattern = "csv", full = T)
listdata_light = lapply(fnames[-7], read.csv, header = T, skip = 1)
light_ID = read.csv(fnames[7], header = T)
names(listdata_light) = paste("", str_extract(basename(fnames[-7]), "Inoue_[a-z]{2,4}_[0-9]{1,6}_[LMH]"), sep = "")
listdata_light = lapply(listdata_light, function(x){x %>% select(Time = 2, Light_lux = 4)})

#データの編集--------------------------------------------------------------------------------------------------
light = listdata_light %>%
  bind_rows(.id = "filename") %>% 
  tidyr::separate(col = filename, into = c("Inoue", "Position", "a", "Density"), sep = "_") %>%
  select(Position, Time, Light_lux, Density) %>% 
  mutate(Time = as.POSIXct(parse_date_time(Time, orders = "%m/%d/%y %I:%M:%S %p", tz = "JST", locale = "ja_JP.utf8")))

light_ID = light_ID %>% 
  mutate(starttime = as.POSIXct(as.character(Starttime), tz = "JST")) %>% 
  select(-3)

light$Hz = NA
light$Density = NA

removetime = 30 #s 初めの切り取り時間
extracttime = removetime + 20 #s 使用する時間
for (i in 1:nrow(light_ID)) {
  
  light[light$Time >= light_ID$starttime[i] + removetime &
          light$Time <= light_ID$starttime[i] + extracttime, "Hz"] = light_ID$Hz[i]
  
  light[light$Time >= light_ID$starttime[i] + removetime &
          light$Time <= light_ID$starttime[i] + extracttime, "Density"] = paste(light_ID$Density[i])
}

light = light %>%
  na.omit() %>% 
  group_by(Hz, Density, Position) %>% 
  mutate(Elapsed_time = as.integer((Time - min(Time)))) %>% 
  ungroup() %>% 
  mutate(Density = recode(Density, H = "High", M = "Middle", L = "Low"))
         

#平均値、ばらつきを計算-------------------------------------------------------
light_mean = light %>% group_by(Density, Hz) %>% 
  summarise(Light_mean = mean(Light_lux), Light_sd = sd(Light_lux))

藻場の形状変化データの編集

## 海藻の面積データ（cm^2）------------------------------------------
fname = dir("../回流水槽光合成実験5/Data/Area/", pattern = "area_cm.csv", full = T)
Area_cm = lapply(fname, read.csv, header = T)

#データの編集
names(Area_cm) = basename(fname)

Area_cm = Area_cm %>%
  bind_rows(.id = "filenames") %>% 
  tidyr::separate(filenames, into = c("Density", "a", "b"), sep = "_") %>% 
  select(Density, Area_cm = "Area") %>% 
  mutate(Hz = rep(c(0,0.5,1,2,4,6,8,10,15), 3),
         Density = recode(Density, H = "High", M = "Middle", L = "Low")) %>% 
  group_by(Density) %>% 
  mutate(Area_cm_rate = (1-(Area_cm / first(Area_cm)))*100) #藻場面積の減少率、単位は％

光合成速度、光、藻場の形状変化、上流流速、藻場の深度グループごとの流速、藻場の高さデータの結合

#藻場（26cm, 41cm）の中（in）、上の淵(edge)、上(above)ごとの流速
Velocity_bed = vector_mean2 %>%
  ungroup() %>% 
  filter(Distance %in% c(41, 71)) %>%
  mutate(Velocity = sqrt((u_mean^2 + v_mean^2 + w_mean^2)/3)) %>%
  group_by(Hz, Density, Depth2, upstream_u_mean) %>%
  summarise(Velocity_bed = mean(Velocity),
            Velocity_bed_sd = sd(Velocity)) 

#藻場（26cm, 41cm）の中（in）、上の淵(edge)、上(above)ごとのTKE
TKE_bed = vector_zansa_mean2 %>%
  ungroup() %>% 
  filter(Distance %in% c(41, 71)) %>%
  group_by(Hz, Density, Depth2, upstream_u_mean) %>%
  summarise(TKE_bed = mean(TKE),
            TKE_bed_sd = sd(TKE)) 

#藻場（26cm, 41cm）の中（in）、上の淵(edge)、上(above)ごとのstress
stress_bed = vector_zansa_mean2 %>%
  ungroup() %>% 
  filter(Distance %in% c(41, 71)) %>%
  group_by(Hz, Density, Depth2, upstream_u_mean) %>%
  summarise(stress_bed = mean(stress),
            stress_bed_sd = sd(stress)) 

all_data = Slope_DO %>% 
  full_join(light_mean, by = c("Hz", "Density")) %>% 
  full_join(Area_cm, by = c("Hz", "Density")) %>% 
  full_join(upstream_u_mean, by = c("Hz", "Density")) %>% 
  full_join(Velocity_bed %>%
              select(-upstream_u_mean) %>% 
              spread(key = Depth2, value = Velocity_bed) %>% 
              rename(Velocity_above = above, Velocity_edge = edge, Velocity_in = 'in'),
            by = c("Hz", "Density")) %>% 
  full_join(stress_bed %>% 
              select(-upstream_u_mean) %>% 
              spread(key = Depth2, value = stress_bed) %>% 
              rename(stress_above = above, stress_edge = edge, stress_in = 'in'),
            by = c("Hz", "Density")) %>% 
  full_join(TKE_bed %>% 
              select(-upstream_u_mean) %>% 
              spread(key = Depth2, value = TKE_bed) %>% 
              rename(TKE_above = above, TKE_edge = edge, TKE_in = 'in'),
            by = c("Hz", "Density"))

#藻場の高さcm
bed_height = read_csv("../回流水槽光合成実験5/Data/Exp_Height_180522.csv") %>% mutate(Height = Height + 1)

all_data = all_data %>% 
  full_join(bed_height, by = c("Hz", "Density")) 

滞留時間と藻場の空隙率（1-solid_volume_fractions）の算出

流速の変化により藻場が変形するため、藻場の体積は変化する。また、密度によって藻場内の海藻の体積も変化する。この2点を滞留時間の計算に用いる。

Length = (90-26)/100 #m 藻場の長さ
Width = 0.30 #m 藻場の幅

all_data = all_data %>% 
  mutate(Total_bed_volume = Height/100 * Length * Width,
         Net_bed_volume = Total_bed_volume - Seaweed_volume_m3,
         Tr_mean = Total_bed_volume / (Velocity_in * Height * Width), #滞留時間の単位はｓ
         Solid_volume_fractions = Net_bed_volume/Total_bed_volume*100)#空隙率単位は％ 

レイノルズ数の算出

レイノルズ数は、流体と物体大きさから、流れの状態を表現する値である。
\(Re=v*L/u\)
\(v: Velocity\)
\(L: 物体の長さ\)
\(u: 動粘度係数\)

#レイノルズ数の計算-------------------------------------
#Re = vL/u v:velocity, L:物体の大きさ
#計算に用いる物体のスケールLは藻場の高さ（m）, 海藻間の距離(m)とする。
u = 0.94*10^-6 #m/s 海水２５度の時の動粘度係数(Bilger and Atkinson 1992, Kays and Crawford 1993, Thomas and Atkinson 1997)

all_data = all_data %>% 
  mutate(Re_height = ((upstream_u_mean) * (Height))/u,
         Distance_individuals = recode(Density, "High" = 3, "Middle" = 6, "Low" = 8),#個体間の距離cm
         Re_distance = ((upstream_u_mean) * (Distance_individuals/100))/u,
         Re_height_distance = (upstream_u_mean * ((Distance_individuals+Height)/100))/u) %>% 
  return()

#密度の順位変更
all_data$Density = factor(all_data$Density, levels = c("High", "Middle", "Low"))

藻場の流動環境の結果

藻場のleading edgeを0cmとし、藻場の終わりは、藻場密度がHigh:54cm, Middle:57cm, Low:54cmとした。

#藻場の高さ（流速、密度条件によって変動する）を追加
vector_mean2 = vector_mean2 %>% 
  left_join(bed_height) %>% 
  mutate(`H/h` = Depth/Height,
         Density = factor(Density, levels = c("High", "Middle", "Low")))

vector_zansa_mean2 = vector_zansa_mean2 %>% 
  left_join(bed_height) %>% 
  mutate(`H/h` = Depth/Height,
         Density = factor(Density, levels = c("High", "Middle", "Low")))

#藻場の縁を0cmとする
vector_mean2 = vector_mean2 %>% 
  mutate(Distance_leading_edge = Distance - 28) #leading_edgeは28cm 

vector_zansa_mean2 = vector_zansa_mean2 %>% 
  mutate(Distance_leading_edge = Distance - 28) #leading_edgeは28cm 

#流速条件（Hz）を上流流速の平均値に置き換えた列を作成
Vel_treatment = upstream_u_mean %>% 
  ungroup() %>% 
  group_by(Hz) %>% 
  summarise(Vel_treatment = round(mean(upstream_u_mean), digits = 2))

vector_mean2 = vector_mean2 %>% 
  left_join(Vel_treatment)

vector_zansa_mean2 = vector_zansa_mean2 %>% 
  left_join(Vel_treatment)

その他の要因

光合成速度、呼吸速度の結果

総光合成速度と他の要因

Stanのコード

functions{
  vector modelcurve(real b0, real b1, real b2, real b3, vector speed) {
    return (b0 + b1 * (1-exp(-b2 / b1 * speed)) .* exp(-b3 / b1 * speed));
  }
}
data {
  int<lower=0> N;                 // データの数
  int<lower=0> M;                 // 期待値に使う結果の数
  vector[N] speed;                // 観測した速度
  vector[N] gross;                  // 観測した総一次生産量
  vector[4] prior_log_sigma;      // 各パラメータの事前分散
  vector[4] prior_log_mu;         // 各パラメータの事前期待値
}
transformed data {
  vector[N] gross10;
  vector[M] speed_predict;
  real speed_increment;
  gross10 = gross*10;
  speed_increment = max(speed) / M;
  speed_predict[1] = 0;
  
  for(m in 2:M) {
    speed_predict[m] = speed_predict[m-1] + speed_increment;
  }
}
parameters {
  real<lower=0> offset;
  real<lower=0> pmax;
  real<lower=0> alpha;
  real<lower=0> beta;
  
  real<lower=0> sigma;
}
transformed parameters {
  vector[N] fitted;
  fitted = modelcurve(offset, pmax, alpha, beta, speed);
}
model {
  offset ~ lognormal(prior_log_mu[1], prior_log_sigma[1]);
  pmax   ~ lognormal(prior_log_mu[2], prior_log_sigma[2]);
  alpha  ~ lognormal(prior_log_mu[3], prior_log_sigma[3]);
  beta   ~ lognormal(prior_log_mu[4], prior_log_sigma[4]);
  sigma ~ cauchy(0, 2.5);
  gross10 ~ lognormal(log(fitted), sigma);
}
generated quantities {
  vector[M] expected10;
  vector[M] expected;
  vector[N] log_lik;
  expected10 = modelcurve(offset, pmax, alpha, beta, speed_predict);
  expected = expected10/10;
  for(n in 1:N) {
    log_lik[n] = lognormal_lpdf(log(gross10[n]) | log(fitted[n]), sigma);
  }
}

モデル式へのあてはめ

#パッケージの読み込み
library(rstan)

#Stanに渡すデータを準備(Stan側と同じオブジェクト名にする)
Gross_data = all_data %>%
  filter(Hz != 0, type == "GP")

N = nrow(Gross_data)                                      #データの数
speed = Gross_data$upstream_u_mean                        #観測した速度
gross = Gross_data$Slope                                  #観測した総一次生産量 

#密度条件ごとに変化しないパラメーター
M = 50                                                  #期待値に使う結果の数
prior_log_sigma = as.vector(c(1,1,1,1))               #各パラメータの事前分散
prior_log_mu = as.vector(c(1,1,1,1))

#Stanの実行
rng_seed　= 1111
nchains = 4
iter = 1000
ncores = 4 

stanout = sampling(nonlinearmodel_gross, 
                   chains = nchains, 
                   iter = iter, 
                   cores = ncores,
                   seed = rng_seed,
                   control = list(adapt_delta = 1))

saveRDS(stanout, file ="../回流水槽光合成実験5/Modified_Data/stanout.obj")　#解析結果が保存されたオブジェクトを保存する

Stanの結果

stanout@model_pars

##  [1] "offset"     "pmax"       "alpha"      "beta"       "sigma"     
##  [6] "fitted"     "expected10" "expected"   "log_lik"    "lp__"

pars =  c( "offset", "pmax", "alpha", "beta", "sigma")
pars_fitted =  c( "fitted")
pars_expected =  c( "expected")
traceplot(stanout, pars = c(pars))

## Stan model '8d896e915aa8a88ef1bd2ae55bf0ec2d' does not contain samples.

# traceplot(stanout, pars = c(pars_fitted))
# traceplot(stanout, pars = c(pars_expected))
print(stanout, prob=c(0.025, 0.5, 0.975), par = pars)

## Stan model '8d896e915aa8a88ef1bd2ae55bf0ec2d' does not contain samples.